rm(list=ls())

library(foreign)
library(tidyr)
library(plyr)
library(sandwich)
library(lmtest)
library(survival)
library(xtable)
library(ggplot2)
library(scales)
library(stargazer)

prop.res <- read.csv("Data/beim_rader_data.csv")

########### plot num Dems and num Reps in prec cases per split

plotdata <- subset(prop.res, select = c(num_dem_prec, num_rep_prec))

pdf("Graphs/num_dem_rep.pdf",height=5,width=5)
ggplot(plotdata, aes(x=jitter(num_dem_prec), y=jitter(num_rep_prec))) +
  geom_point(size=2, shape=1) + 
  geom_hline(yintercept = median(plotdata$num_rep_prec), lty="dashed") +
  geom_vline(xintercept = median(plotdata$num_dem_prec), lty="dashed") +
  theme_bw() + xlab("Number of Democratic Appointees") + ylab("Number of Republican Appointees") +
  scale_x_continuous(breaks=pretty_breaks()) + scale_y_continuous(breaks=pretty_breaks()) + 
  theme(plot.title = element_text(size = 13, hjust=.5),
        plot.margin =       unit(c(1.3, 1.3, 1, 1), "cm"),
        axis.text.y = element_text(size=11), axis.text.x = element_text(size=11), 
        axis.title.x = element_text(size=10, vjust=-1), axis.title.y = element_text(size=10, vjust=3))
dev.off()

############## Polarization histogram
plotdata <- subset(prop.res, select=c(polarization))

pdf("Graphs/polarization_plot.pdf",height=5,width=5.5)
ggplot(plotdata, aes(polarization)) + 
  geom_histogram(color="black",fill="gray70",breaks=seq(0, 1, by = .1)) +
  geom_vline(xintercept = mean(plotdata$polarization,na.rm=T),color="blue", lwd=1) +
  geom_vline(xintercept = median(plotdata$polarization,na.rm=T), color="blue", lwd=1, lty="dashed") +
  theme_bw() + xlab("Polarization") + ylab("Splits") +
  scale_x_continuous(breaks=pretty_breaks()) +  
  theme(plot.title = element_text(size = 13, hjust=.5),
        plot.margin =       unit(c(1.3, 1.3, 1, 1), "cm"),
        axis.text.y = element_text(size=11), axis.text.x = element_text(size=11), 
        axis.title.x = element_text(size=12, vjust=-1), axis.title.y = element_text(size=12, vjust=3))
dev.off()

################ Num circ histogram
pdf("Graphs/numcirc_density_by_resolved.pdf")
plot(density(prop.res$num_circ[prop.res$resolved==1]), ylim=c(0,.3), xlab="Number of Circuits", main="")
lines(density(prop.res$num_circ[prop.res$resolved==0]), col="red")
text(2.5,.19, "Unresolved", col="red")
text(10,.09, "Resolved")
abline(v=median(prop.res$num_circ[prop.res$resolved==1]))
abline(v=median(prop.res$num_circ[prop.res$resolved==0]), col="red")
dev.off()

##############
############## MAIN REGRESSION
############## split level polarization LPM, with robust SE (HC3 is default type)
##############

lpm.pol <- lm(resolved~polarization + z_num_circ + dormant_start + team_size_circ_diff +
                any_dissent  + any_sgpet + crim + econ, data=prop.res)
lpm.pol$se<-vcovHC(lpm.pol)

########## COEF PLOT, split level polarization LPM #########
plotdata<-cbind.data.frame(coeftest(lpm.pol, vcov=lpm.pol$se)[-1,1],
                           coeftest(lpm.pol, vcov=lpm.pol$se)[-1,2],
                           c("Polarization", "Number of Circuits","Dormant Start","Lopsidedness",
                             "Any Dissents", "Any SG petitions", "Criminal Procedure","Economic Activity"))
   colnames(plotdata) <- c("Coef","SE","Var")
plotdata$Var <- factor(plotdata$Var, levels = c("Polarization", "Number of Circuits","Dormant Start","Lopsidedness",
                                                "Any Dissents", "Any SG petitions", "Criminal Procedure","Economic Activity"))

pdf("Graphs/coefplot_lpm_polarization.pdf",height=5,width=5)
ggplot(plotdata, aes(x=Var, y=Coef, ymin = Coef - qnorm(1 - 0.05 / 2) * SE,
                     ymax = Coef + qnorm(1 - 0.05 / 2) * SE)) +
  scale_x_discrete(limits = rev(levels(plotdata$Var))) +
  geom_pointrange(size=.55) +      
  geom_hline(yintercept = 0, lwd = .75, lty="dotted", alpha = .4) +
  coord_flip() + theme_bw() +
  theme(plot.title = element_text(size = 13, hjust=.5),
        strip.text = element_text(size=11, vjust=.7), panel.spacing = unit(.75, "lines"),
        plot.margin =       unit(c(1, 1, 0.5, 0.5), "lines"),
        strip.background = element_rect(colour="white", fill="white"), 
        axis.text.y = element_text(size=12), axis.text.x = element_text(size=11), 
        axis.title.x = element_blank(), axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.major =  element_line(colour = "grey90", size = 0.2))
dev.off()


########## TABLE ###########
lpm_pol_tab <- stargazer(lpm.pol, type = "latex", digits=2,
                         se = list(sqrt(diag(lpm.pol$se))),
                         covariate.labels = c("Split Polarization", "Number of Circuits","Dormant Start","Lopsidedness",
                                              "Any Dissents", "Any SG petitions", "Criminal Procedure","Economic Activity","Intercept"),
                         model.numbers = FALSE, 
                         dep.var.labels = c("Split Resolved"),
                         star.char = c("", "*","**"), 
                         star.cutoffs = c(0.2, 0.1, 0.05), 
                         label = "table:pol_lpm")
note.latex <- "\\multicolumn{2}{l} {\\parbox[t]{.65\\textwidth}{ \\textit{Notes:} $^{*}$p$<$0.1; $^{**}$p$<$0.05; 
Number of Circuits is standardized. HC3 robust standard errors shown.}} \\\\"

lpm_pol_tab[grepl("Note",lpm_pol_tab)] <- note.latex
cat(lpm_pol_tab, sep = "\n")
write(lpm_pol_tab, "Graphs/lpm_pol.txt")



##############
############## supplementary regressions
############## with robust SE (HC3 is default type)
##############

##### SG pet and polarization

lpm.sg <- lm(any_sgpet~polarization, data=prop.res)
lpm.sg$se<-vcovHC(lpm.sg)


lpm.sg2 <- lm(any_sgpet~polarization + z_num_circ + dormant_start + team_size_circ_diff +
               any_dissent + crim + econ, data=prop.res)
lpm.sg2$se<-vcovHC(lpm.sg2)

sg_pol_tab <- stargazer(lpm.sg, lpm.sg2, type = "latex", digits=2,
                        label = "table:sg_pol_lpm",
                        se = list(sqrt(diag(lpm.sg$se)), sqrt(diag(lpm.sg2$se))),
                        covariate.labels = c("Split Polarization", "Number of Circuits","Dormant Start","Lopsidedness",
                                              "Any Dissents", "Criminal Procedure","Economic Activity","Intercept"),
                         model.numbers = FALSE, 
                         dep.var.labels = c("SG Petition"),
                         star.char = c("", "*","**"), 
                         star.cutoffs = c(0.2, 0.1, 0.05))
note.latex <- "\\multicolumn{3}{l} {\\parbox[t]{.75\\textwidth}{ \\textit{Notes:} $^{*}$p$<$0.1; $^{**}$p$<$0.05; 
Number of Circuits is standardized. HC3 robust standard errors shown.}} \\\\"

sg_pol_tab[grepl("Note",sg_pol_tab)] <- note.latex
cat(sg_pol_tab, sep = "\n")
write(sg_pol_tab, "Graphs/sg_pol.txt")


###Granted, rather than resolved
lpm.pol.grant <- lm(petgrant~polarization + z_num_circ + dormant_start + team_size_circ_diff +
                any_dissent  + any_sgpet + crim + econ, data=prop.res)
lpm.pol.grant$se<-vcovHC(lpm.pol.grant)

granted_table <- stargazer(lpm.pol.grant, type="latex", digits=2,
							label= "table:granted_lpm",
							se = list(sqrt(diag(lpm.pol.grant$se))),
							covariate.labels = c("Split Polarization", "Number of Circuits", "Dormant Start", "Lopsidedness",
							"Any Dissents", "Criminal Procedure", "Economic Activity", "Intercept"),
							model.numbers= FALSE,
							dep.var.labels= c("Certiorari Granted"),
							star.char=c("", "*", "**"),
							star.cutoffs = c(0.2, 0.1, 0.05))
note.latex <- "\\multicolumn{3}{l} {\\parbox[t]{.75\\textwidth}{ \\textit{Notes:} $^{*}$p$<$0.1; $^{**}$p$<$0.05; 
Number of Circuits is standardized. HC3 robust standard errors shown.}} \\\\"
lpm.pol.grant[grepl("Note",granted_table)] <- note.latex
cat(granted_table, sep = "\n")
write(granted_table, "Graphs/grantedtable.txt")




##############
############## LOGIT VERSION FOR APPENDIX
############## split level polarization LPM, with robust SE (HC3 is default type)
##############

log.pol <- glm(resolved~polarization + z_num_circ + dormant_start + team_size_circ_diff +
                any_dissent  + any_sgpet + crim + econ, data=prop.res, family = "binomial")

########## COEF PLOT, split level polarization logit #########
plotdata<-cbind.data.frame(summary(log.pol)$coefficients[-1,1],
                           summary(log.pol)$coefficients[-1,2],
                           c("Polarization", "Number of Circuits","Dormant Start","Lopsidedness",
                             "Any Dissents", "Any SG petitions", "Criminal Procedure","Economic Activity"))
colnames(plotdata) <- c("Coef","SE","Var")
plotdata$Var <- factor(plotdata$Var, levels = c("Polarization", "Number of Circuits","Dormant Start","Lopsidedness",
                                                "Any Dissents", "Any SG petitions", "Criminal Procedure","Economic Activity"))

pdf("Graphs/coefplot_logit_polarization.pdf",height=5,width=5)
ggplot(plotdata, aes(x=Var, y=Coef, ymin = Coef - qnorm(1 - 0.05 / 2) * SE,
                     ymax = Coef + qnorm(1 - 0.05 / 2) * SE)) +
  scale_x_discrete(limits = rev(levels(plotdata$Var))) +
  geom_pointrange(size=.55) +      
  geom_hline(yintercept = 0, lwd = .75, lty="dotted", alpha = .4) +
  coord_flip() + theme_bw() +
  theme(plot.title = element_text(size = 13, hjust=.5),
        strip.text = element_text(size=11, vjust=.7), panel.spacing = unit(.75, "lines"),
        plot.margin =       unit(c(1, 1, 0.5, 0.5), "lines"),
        strip.background = element_rect(colour="white", fill="white"), 
        axis.text.y = element_text(size=12), axis.text.x = element_text(size=11), 
        axis.title.x = element_blank(), axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.major =  element_line(colour = "grey90", size = 0.2))
dev.off()

########## TABLE LOGIT FOR APPENDIX ###########

log_pol_tab <- stargazer(log.pol, type = "latex", digits=2,
                          label = "table:pol_logit",
                          covariate.labels = c("Split Polarization", "Number of Circuits","Dormant Start","Lopsidedness",
                                               "Any Dissents", "Any SG petitions","Criminal Procedure","Economic Activity","Intercept"),
                          model.numbers = FALSE, 
                          dep.var.labels = c("Split Resolved"),
                          star.char = c("", "*","**"), 
                          star.cutoffs = c(0.2, 0.1, 0.05))
note.latex <- "\\multicolumn{2}{l} {\\parbox[t]{.75\\textwidth}{ \\textit{Notes:} $^{*}$p$<$0.1; $^{**}$p$<$0.05; 
Number of Circuits is standardized.}} \\\\"
log_pol_tab[grepl("Note",log_pol_tab)] <- note.latex
cat(log_pol_tab, sep = "\n")
write(log_pol_tab, "Graphs/log_pol.txt")


############################ 
############################  
############## Does split polarization predict SC polarization?
############################ 
############################ 

prop.res$sc_unan <- ifelse(prop.res$connected=="unan",1,0)

### regressions using SCpolarization, which is analogous to polarization
# unan = 0
SCpolarization <- lm(SCpolarization~polarization + z_num_circ + dormant_start + team_size_circ_diff +
             any_dissent + any_sgpet + crim + econ, data=prop.res, subset=resolved==1)

### regressions using connected
prop.res$connected0 <- ifelse(prop.res$connected=="connected",1,0)
# unan = 0
SCconnected <- lm(connected0~polarization + z_num_circ + dormant_start + team_size_circ_diff +
             any_dissent + any_sgpet + crim + econ, data=prop.res, subset=resolved==1)

### regressions using sc_SCpolarization, diff in mean SC scores
prop.res$sc_SCpolarization0 <- ifelse(is.na(prop.res$sc_SCpolarization),0,prop.res$sc_SCpolarization)
# unan = 0
SCSC <-lm(sc_SCpolarization0~polarization + z_num_circ + dormant_start + team_size_circ_diff +
             any_dissent + any_sgpet + crim + econ, data=prop.res, subset=resolved==1)

### regressions using sc_MQpolarization, diff in mean MQ scores
prop.res$mq_SCpolarization0 <- ifelse(is.na(prop.res$mq_SCpolarization),0,prop.res$mq_SCpolarization)
prop.res$z_mq_SCpolarization0 <- (prop.res$mq_SCpolarization0 - mean(prop.res$mq_SCpolarization0))/(2*sd(prop.res$mq_SCpolarization0))


# unan = 0
SCMQ <- lm(mq_SCpolarization0~polarization + z_num_circ + dormant_start + team_size_circ_diff +
                     any_dissent + any_sgpet + crim + econ, data=prop.res, subset=resolved==1)


###################### SC polarization Table with three regressions
###################### 



### FYI make sure to add to latex document preamble \usepackage[usestackEOL]{stackengine}

sc_pol_tab <- stargazer(SCpolarization, SCMQ, SCconnected, type = "latex", digits=2,
                        label = "table:sc_polarization",
                        covariate.labels = c("Split Polarization", "Number of Circuits","Dormant Start","Lopsidedness",
                               "Any Dissents", "Any SG petitions", "Criminal Procedure","Economic Activity","Intercept"),
              column.labels=c("\\addstackgap[5pt]{\\shortstack{Analogous polarization\\\\ measure}}", 
                              "\\shortstack{Difference in mean\\\\  Martin-Quinn score}", 
                              "\\shortstack{Connected\\\\  coalitions}"), 
              model.numbers = FALSE, dep.var.labels.include = FALSE,
              star.char = c("", "*","**"),
              star.cutoffs = c(0.2, 0.1, 0.05))
note.latex <- "\\multicolumn{4}{l} {\\parbox[t]{\\textwidth}{ \\textit{Notes:} $^{*}$p$<$0.1; $^{**}$p$<$0.05; 
Number of Circuits is standardized.}} \\\\"

sc_pol_tab[grepl("Note",sc_pol_tab)] <- note.latex
cat (sc_pol_tab, sep = "\n")
write(sc_pol_tab, "Graphs/sc_polarization.txt")


###################### SC polarization coef graph with three regressions
###################### 

plotdata<-cbind.data.frame(c(summary(SCpolarization)$coefficients[2,1],summary(SCMQ)$coefficients[2,1],summary(SCconnected)$coefficients[2,1]),
                           c(summary(SCpolarization)$coefficients[2,2],summary(SCMQ)$coefficients[2,2],summary(SCconnected)$coefficients[2,2]),
                           c(rep("Polarization",3)),
                           c("Analogous polarization\nmeasure","Difference in mean\nMartin Quinn score","Connected\ncoalitions"))
colnames(plotdata) <- c("Coef","SE","Var","DV")
plotdata$DV <- factor(plotdata$DV, levels = c("Analogous polarization\nmeasure","Difference in mean\nMartin Quinn score","Connected\ncoalitions"))

pdf("Graphs/coefplot_sc_polarization.pdf",height=4,width=2.5)
ggplot(plotdata, aes(x=Var, y=Coef, ymin = Coef - qnorm(1 - 0.05 / 2) * SE,
                     ymax = Coef + qnorm(1 - 0.05 / 2) * SE)) +
  scale_x_discrete(limits = rev(levels(plotdata$Var))) +
  xlab("Polarization\n") +
  facet_wrap(~ DV, ncol=1) +
  #facet_wrap(~ DV) +
  geom_pointrange(size=.55) +      
  geom_hline(yintercept = 0, lwd = .75, lty="dotted", alpha = .4) +
  coord_flip() + theme_bw() +
  theme(strip.text = element_text(size=8, vjust=-1), panel.spacing = unit(.2, "lines"),
        plot.margin =       unit(c(1, 1, 0.5, 1), "lines"),
        #plot.margin =       unit(c(1, 1, 1, 3), "lines"),
        strip.background = element_rect(colour="white", fill="white"), 
        axis.text.y = element_blank(), axis.text.x = element_text(size=11), 
        axis.title.x = element_blank(), axis.title.y = element_text(size=13),
        axis.ticks.y = element_blank(),
        panel.grid.major =  element_line(colour = "grey90", size = 0.2))
dev.off()